/*!
 ***************************************************************************
 * \file heat1.cpp
 * Solving the Heat difusion equation in 1D.
 ***************************************************************************
 *  \author    Luis M. de la Cruz Salas [ Mon Mar 17 08:47:37 GMT 2008 ]
 ***************************************************************************
 * This code solves the unsteady heat diffusion equation in one dimension:
 * \f$ \frac{\partial T}{\partial t} = \nabla^2 T \f$.
 * This equation is solved in the interval \f$ x \in [0,1] \f$.
 * The initial condition for the temperature is \f$ T = 10 \cos (3 \pi x)\f$,
 * and the boundary conditions:  \f$ T=10 \f$ for \f$ x=0 \f$, and \f$ T=-10 \f$ 
 * for \f$ x = 1 \f$.
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 *  \verbatim 
   % make
   % ./heat1  \endverbatim
 * The results are stored in \c Data1
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"
#include "Schemes/CDS.hpp"
#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,1>::huge ScalarField1D;

int main() 
{
    double length, dt, tolerance; 
    int num_nodes, frequency;
    double Amplitude;

// Reading data from the input file
    std::ifstream input_cin ("input1"); 
    input_cin >> length
	      >> num_nodes
	      >> Amplitude
	      >> dt
	      >> tolerance
	      >> frequency;
    input_cin.close();

    StructuredMesh<Uniform<double, 1> > mesh(length, num_nodes);
    float dx = mesh.getDelta(X);
    mesh.print();
    ScalarField1D T ( mesh.getExtentVolumes() ); 
    ScalarField1D Tn( mesh.getExtentNodes()   );

    firstIndex i;
    Tn = Amplitude * cos (3 * PI * dx * i);
    NumUtils::interpolateToVolumes(T, Tn);

#ifdef WITH_DX
    InOut::writeToFile_DX(Tn, 0, "./Data1/temp.", dx);
#endif

#ifdef WITH_GNUPLOT
    InOut::writeToFile(Tn, 0, "./Data1/temp.", dx);
    GNUplot plotter;
    plotter("set grid");
    plotter("p \"./Data1/temp.00000\" title \"Initial cond.\"w lp");
    std::cout << "\n >---> Press <1> and then <enter> to continue \n";
    int pause; std::cin >> pause;
#endif

    SparseMatrix< Diagonal<double, 1> > A(num_nodes); 

    ScalarField1D              b(num_nodes);

    ScalarEquation< CDS<double, 1> > heat(T, A, b, mesh.getDeltas());
    heat.setDeltaTime(dt); 
    heat.setGamma(1.0);
    heat.setDirichlet(LEFT_WALL); 
    heat.setDirichlet(RIGHT_WALL);
    heat.print();

    ofstream error_fp("./Data1/error"), residual_fp("./Data1/residual");
    int iteration = 0;
    double error = 1.0;
    double residual = 1.0;

    while (error > tolerance) {
      std::cout << "\n ---- Iteration " << ++iteration 
		<< "\t Res = " << residual
		<< "\t Err = " << error;
      heat.calcCoefficients();
      Solver::TDMA1D(heat);
      residual = heat.calcResidual();
      error = NumUtils::calcErrorL1(T, heat.phi);

      error_fp << iteration << "\t" << error << std::endl;
      residual_fp << iteration << "\t" << residual << std::endl;	
 
      heat.update();
      if ( !(iteration % frequency) ) {
	NumUtils::interpolateToNodes(Tn, T);

#ifdef WITH_DX
	InOut::writeToFile_DX(Tn, iteration, "./Data1/temp.", dx);
#endif
	
#ifdef WITH_GNUPLOT
	InOut::writeToFile(Tn, 1, "./Data1/temp.",dx);
	std::cout << "\n >---> Press <1> and then <enter> to continue \n";
	std::cin >> pause;	
	plotter("p \"./Data1/temp.00000\" title \"Initial cond.\" w lp, \"Data1/temp.00001\" title \"Num. sol.\" w lp");
#endif
      }
    }

#ifdef WITH_DX
	std::cout << "\n\n NOTE: The data were written in the [./Data1] directory in DX format";
#endif

    return 0;
}




